Load phyloseq data generated from dada2

# Load Phyloseq object generated in HIV_oral_preprocessing.Rmd
# Options include GreenGenes, SILVA, RDP and HitDB
ps0 <- readRDS("~/Dropbox/Research/human_volunteer/data/ps0.Harris.rdp.RDS")
ps0
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1936 taxa and 122 samples ]
## tax_table()   Taxonomy Table:    [ 1936 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1936 tips and 1934 internal nodes ]
# Load mapping file
map <- import_qiime_sample_data("~/Dropbox/Research/human_volunteer/data/mapping_human_volunteer.txt")
ps0 <- merge_phyloseq(ps0, map)
ps0
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1936 taxa and 122 samples ]
## sample_data() Sample Data:       [ 122 samples by 119 sample variables ]
## tax_table()   Taxonomy Table:    [ 1936 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1936 tips and 1934 internal nodes ]
# View sample variables & generate basic stats
sample_variables(ps0)
##   [1] "X.SampleID"                  "BarcodeSequence"            
##   [3] "LinkerPrimerSequence"        "sequence_ID"                
##   [5] "miseq_run"                   "plate"                      
##   [7] "well"                        "collaborator"               
##   [9] "tissue"                      "number_seqs_slout"          
##  [11] "sample_surviving_dada2"      "sample_ID"                  
##  [13] "tube_ID"                     "patient_ID"                 
##  [15] "randomization"               "day"                        
##  [17] "randomization_arm"           "Date_assessment"            
##  [19] "Age"                         "Gender"                     
##  [21] "Ethnicity"                   "Asthma"                     
##  [23] "Eczema"                      "Hayfever"                   
##  [25] "Allergies"                   "Def_allergies"              
##  [27] "Other_illness"               "Chronic_medication"         
##  [29] "Def_medication"              "Herbs_OTC"                  
##  [31] "def_Herbs_OTC"               "Alcohol"                    
##  [33] "Alcoholunitsweek"            "Smoking"                    
##  [35] "Drugs"                       "Drugstype"                  
##  [37] "Drugsusage"                  "Weight"                     
##  [39] "Length"                      "BMIcalc"                    
##  [41] "SystolicBP"                  "DiastolicBP"                
##  [43] "Pulse"                       "Respfreq"                   
##  [45] "Physical_exam_abnorm"        "def_phys_abnormalities"     
##  [47] "Lab_abnorm_screen"           "def_abnorm_lab_screen"      
##  [49] "Serum_screen"                "PBMC_screen"                
##  [51] "Feces_screen"                "Date_day_9"                 
##  [53] "New_medication_day_9"        "def_new_med_day_9"          
##  [55] "Course_completed"            "Antibiotic_side_effects"    
##  [57] "Def_abx_side_effects"        "Lab_abnorm_post_antibiotics"
##  [59] "def_abnorm_post_antibiotics" "Date_day_0"                 
##  [61] "Feces_day_0"                 "New_medication_day0"        
##  [63] "Def_new_medication_day_0"    "Vaccination_side_effects"   
##  [65] "Def_vac_side_effects"        "Lab_abnorm_post_vaccination"
##  [67] "def_abnorm_post_vaccination" "Date_day_7"                 
##  [69] "Feces_day_1"                 "Feces_day_2"                
##  [71] "Feces_day_3"                 "Feces_day_4"                
##  [73] "Feces_day_5"                 "Feces_day_6"                
##  [75] "Feces_day_7"                 "Serum_day_7"                
##  [77] "New_medication_day7"         "def_new_med_day7"           
##  [79] "Date_day_14"                 "New_medication_day14"       
##  [81] "def_new_med_day14"           "Serum_day_14"               
##  [83] "PBMC_day_14"                 "Date_day_28"                
##  [85] "New_medication_day28"        "def_new_med_day28"          
##  [87] "Serum_day_28"                "PBMC_day_28"                
##  [89] "Study_complete"              "pre_RV_IgA"                 
##  [91] "d7_RV_IgA"                   "boost_change"               
##  [93] "d7_RV_IgA_boost_fc"          "d7_roto_boost"              
##  [95] "d7_rota_boost_updated"       "d14_RV_IgA"                 
##  [97] "d28_RV_IgA"                  "RV_Ag_shedding_d0"          
##  [99] "RV_Ag_shedding_d1"           "RV_Ag_shedding_d2"          
## [101] "RV_Ag_shedding_d3"           "RV_Ag_shedding_d4"          
## [103] "RV_Ag_shedding_d5"           "RV_Ag_shedding_d6"          
## [105] "Shedding"                    "pneumo_IgG_pre"             
## [107] "pneumo_IgG_7"                "pneumo_IgG_14"              
## [109] "pneumo_IgG_28"               "pneumo_IgG_7vPre"           
## [111] "pneumo_IgG_14vPre"           "pneumo_IgG_28vPre"          
## [113] "tetanus_IgG_pre"             "tetanus_IgG_d7"             
## [115] "tetanus_IgG_d14"             "tetanus_IgG_d28"            
## [117] "tetanus_IgG_7vPre"           "tetanus_IgG_14vPre"         
## [119] "tetanus_IgG_28vPre"
sd(sample_sums(ps0))
## [1] 23974.15
get_taxa_unique(ps0, "Phylum")
##  [1] "Bacteroidetes"             "Firmicutes"               
##  [3] "Verrucomicrobia"           "Proteobacteria"           
##  [5] "Fusobacteria"              "Synergistetes"            
##  [7] "Actinobacteria"            "Spirochaetes"             
##  [9] NA                          "Lentisphaerae"            
## [11] "Euryarchaeota"             "Elusimicrobia"            
## [13] "Cyanobacteria/Chloroplast" "Tenericutes"

ASV summary statistics

# Create a new data frame of the sorted row sums, a column of sorted values from 1 to the total number of individuals/counts for each RSV and a categorical variable stating these are all RSVs.
readsumsdf = data.frame(nreads = sort(taxa_sums(ps0), TRUE), 
                        sorted = 1:ntaxa(ps0),
                        type = "ASVs")


# Add a column of sample sums (total number of individuals per sample)
readsumsdf = rbind(readsumsdf,
                   data.frame(nreads = sort(sample_sums(ps0), TRUE),
                              sorted = 1:nsamples(ps0),
                              type = "Samples"))

# Make a data frame with a column for the read counts of each sample for histogram production
sample_sum_df <- data.frame(sum = sample_sums(ps0))

# Make plots
# Generates a bar plot with # of reads (y-axis) for each taxa. Sorted from most to least abundant
# Generates a second bar plot with # of reads (y-axis) per sample. Sorted from most to least
p.reads = ggplot(readsumsdf, aes(x = sorted, y = nreads)) +
  geom_bar(stat = "identity") +
  ggtitle("ASV Assessment") +
  scale_y_log10() +
  facet_wrap(~type, scales = "free") +
  ylab("# of Reads")

# Histogram of the number of Samples (y-axis) at various read depths
p.reads.hist <- ggplot(sample_sum_df, aes(x = sum)) + 
  geom_histogram(color = "black", fill = "firebrick3", binwidth = 150) +
  ggtitle("Distribution of sample sequencing depth") + 
  xlab("Counts") +
  ylab("# of Samples")

# Final plot, side-by-side
grid.arrange(p.reads, p.reads.hist, ncol = 2)

# Basic summary statistics
summary(sample_sums(ps0))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0    1205   20670   24463   41152   99585

Interpretation: Lots (20+) samples with zero ASV’s. These will likely need to be removed in the subsequent step. Other than that the distributions look normal. Skew from high to low ASV is due to the incorporation of antibiotics treated samples which result in fewer ASV.

Bad sample identification and removal

# Format a data table to combine sample summary data with sample variable data
ss <- sample_sums(ps0)
sd <- as.data.frame(sample_data(ps0)) # useful to coerce the phyloseq object into an R data frame
ss.df <- merge(sd, data.frame("ASV" = ss), by ="row.names") # merge ss with sd by row names. Rename ss to ASVs in the new data frame

# Plot the data by the treatment variable
y = 100 # Set a threshold for the minimum number of acceptable reads. Can start as a guess
x = "randomization_arm" # Set the x-axis variable you want to examine
label = "Description" # This is the label you want to overlay on the points that are below threshold y. Should be something sample specific

# Plot
p.ss.boxplot<- ggplot(ss.df, aes_string(x, y = "ASV")) + # x is what you assigned it above
  geom_boxplot(outlier.colour="NA", aes(group = randomization_arm)) +
  scale_y_log10() +
  geom_hline(yintercept = y, lty = 2) + # Draws a dashed line across the threshold you set above as y
  geom_jitter(alpha = 0.6, width = 0.15, size = 2.5) +
  #geom_text(data = subset(ss.df, RSVs <= y), aes_string(x, y="RSVs", label=label), size=2, hjust -1) +# This labels a subset that fall below threshold variable y and labels them with the label variable
  ggtitle("Number of ASV") +
  ylab("ASV (log10)") +
  facet_grid(~day) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(legend.position = "NULL") +
  theme(axis.title.x = element_blank())
p.ss.boxplot

# Remove samples with fewer than 100 ASV
# Save samples that fell below the selected threshold
ps.failed <- prune_samples(sample_sums(ps0) < y, ps0)
nsamples(ps.failed)
## [1] 22
# Create a table of failed sample information
failed.samples <- data.frame(sample_names(ps.failed), sample_data(ps.failed)$sample_ID, sample_data(ps.failed)$tube_ID)
write.table(failed.samples, file = "./results/failed_samples.txt", sep = "\t")

# Remove samples with fewer than y ASV
ps0
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1936 taxa and 122 samples ]
## sample_data() Sample Data:       [ 122 samples by 119 sample variables ]
## tax_table()   Taxonomy Table:    [ 1936 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1936 tips and 1934 internal nodes ]
ps0 <- prune_samples(sample_sums(ps0) >=y, ps0)
ps0
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1936 taxa and 100 samples ]
## sample_data() Sample Data:       [ 100 samples by 119 sample variables ]
## tax_table()   Taxonomy Table:    [ 1936 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1936 tips and 1934 internal nodes ]
min(sample_sums(ps0))
## [1] 103

Interpretation: There does not appear to be a relationship between low numbers of ASVs, treatment or time-point. So removing samples that behave unexpectedly (< 100 sequences).

Overall sample relationship to evaluate sample outliers

# Outlier evaluation
out.bray <- ordinate(ps0, method = "NMDS", distance = "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1817636 
## Run 1 stress 0.1867117 
## Run 2 stress 0.1883909 
## Run 3 stress 0.1836108 
## Run 4 stress 0.1833127 
## Run 5 stress 0.1889399 
## Run 6 stress 0.1870656 
## Run 7 stress 0.1925691 
## Run 8 stress 0.1880074 
## Run 9 stress 0.1927114 
## Run 10 stress 0.1895174 
## Run 11 stress 0.188155 
## Run 12 stress 0.1821197 
## ... Procrustes: rmse 0.02361735  max resid 0.1442131 
## Run 13 stress 0.1965863 
## Run 14 stress 0.18175 
## ... New best solution
## ... Procrustes: rmse 0.000991212  max resid 0.009441048 
## ... Similar to previous best
## Run 15 stress 0.1919634 
## Run 16 stress 0.1966386 
## Run 17 stress 0.184672 
## Run 18 stress 0.1870022 
## Run 19 stress 0.2003686 
## Run 20 stress 0.1821406 
## ... Procrustes: rmse 0.02418486  max resid 0.1442341 
## *** Solution reached
p.NMDS.outlier <- plot_ordination(ps0, out.bray, color = "randomization_arm", axes = c(1,2)) +
  theme_bw() +
  facet_grid(~day) +
  geom_point(size = 2) +
  ggtitle("NMDS of Bray Distances \nOutlier Evaluation") +
  stat_ellipse(type = "norm", geom = "polygon", alpha = 1/10, aes(fill = randomization_arm))
p.NMDS.outlier

Interpretation: There is no distiguishable feature to identify samples as obvious outliers using this projection. No samples will be removed based on NMDS of Bray-curtis distances.

Prevalence assessment and taxon filtering

# Begin by removing sequences that were not classified as Bacteria or were classified as either mitochondria or chlorplast

ps0 # Check the number of taxa prior to removal
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1936 taxa and 100 samples ]
## sample_data() Sample Data:       [ 100 samples by 119 sample variables ]
## tax_table()   Taxonomy Table:    [ 1936 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1936 tips and 1934 internal nodes ]
ps1 <- ps0 %>%
  subset_taxa(
    Family  != "mitochondria" &
    Class   != "Chloroplast"
  )
ps1 # Confirm that the taxa were removed
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1590 taxa and 100 samples ]
## sample_data() Sample Data:       [ 100 samples by 119 sample variables ]
## tax_table()   Taxonomy Table:    [ 1590 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1590 tips and 1588 internal nodes ]
saveRDS(ps1, file = "./results/ps1.RDS")

# Intermediate data subsetting for prevelance plotting
# Full subsetting with removal of zero count data will be performed in a subsequent chunk
ps1.d_9 <- subset_samples(ps1, day == "-9")
ps1.d0 <- subset_samples(ps1, day == "0")
ps1.d7 <- subset_samples(ps1, day == "7")

# Prevelance estimation
# Calculate feature prevelance across the data set
prevdf.d_9 <- apply(X = otu_table(ps1.d_9),MARGIN = ifelse(taxa_are_rows(ps1.d_9), yes = 1, no = 2),FUN = function(x){sum(x > 0)})
prevdf.d0 <- apply(X = otu_table(ps1.d0),MARGIN = ifelse(taxa_are_rows(ps1.d0), yes = 1, no = 2),FUN = function(x){sum(x > 0)})
prevdf.d7 <- apply(X = otu_table(ps1.d7),MARGIN = ifelse(taxa_are_rows(ps1.d7), yes = 1, no = 2),FUN = function(x){sum(x > 0)})

# Add taxonomy and total read counts to prevdf
prevdf.d_9 <- data.frame(Prevalence = prevdf.d_9, TotalAbundance = taxa_sums(ps1.d_9), tax_table(ps1.d_9))
prevdf.d0 <- data.frame(Prevalence = prevdf.d0, TotalAbundance = taxa_sums(ps1.d0), tax_table(ps1.d0))
prevdf.d7 <- data.frame(Prevalence = prevdf.d7, TotalAbundance = taxa_sums(ps1.d7), tax_table(ps1.d7))

# Create a table of Phylum, their mean abundances across all samples, and the number of samples they were detected in
plyr::ddply(prevdf.d_9, "Phylum", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))})
##             Phylum        1    2
## 1   Actinobacteria 2.000000   92
## 2    Bacteroidetes 2.308594  591
## 3    Elusimicrobia 0.500000    1
## 4    Euryarchaeota 3.000000    6
## 5       Firmicutes 2.067068 2404
## 6     Fusobacteria 0.200000    3
## 7    Lentisphaerae 1.400000    7
## 8   Proteobacteria 1.282609  118
## 9     Spirochaetes 0.000000    0
## 10   Synergistetes 0.000000    0
## 11     Tenericutes 0.000000    0
## 12 Verrucomicrobia 2.800000   14
plyr::ddply(prevdf.d0, "Phylum", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))})
##             Phylum         1    2
## 1   Actinobacteria 1.8478261   85
## 2    Bacteroidetes 1.2265625  314
## 3    Elusimicrobia 0.0000000    0
## 4    Euryarchaeota 3.0000000    6
## 5       Firmicutes 1.1160791 1298
## 6     Fusobacteria 0.4666667    7
## 7    Lentisphaerae 0.8000000    4
## 8   Proteobacteria 2.0543478  189
## 9     Spirochaetes 1.0000000    1
## 10   Synergistetes 1.5000000    3
## 11     Tenericutes 1.0000000    1
## 12 Verrucomicrobia 3.2000000   16
plyr::ddply(prevdf.d7, "Phylum", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))})
##             Phylum        1    2
## 1   Actinobacteria 2.347826  108
## 2    Bacteroidetes 2.007812  514
## 3    Elusimicrobia 0.500000    1
## 4    Euryarchaeota 3.000000    6
## 5       Firmicutes 2.480653 2885
## 6     Fusobacteria 1.133333   17
## 7    Lentisphaerae 1.400000    7
## 8   Proteobacteria 2.239130  206
## 9     Spirochaetes 1.000000    1
## 10   Synergistetes 1.500000    3
## 11     Tenericutes 0.000000    0
## 12 Verrucomicrobia 4.400000   22
# Reorder based on visual inspection of relative dominance of each taxa for plotting
prevdf.d_9$Phylum <- factor(prevdf.d_9$Phylum, levels = c("Firmicutes", "Bacteroidetes", "Proteobacteria", "Actinobacteria", "Verrucomicrobia", "Lentisphaerae", "Fusobacteria", "Elusimicrobia", "Spirochaetes", "Synergistetes", "Tenericutes"))
prevdf.d0$Phylum <- factor(prevdf.d0$Phylum, levels = c("Firmicutes", "Bacteroidetes", "Proteobacteria", "Actinobacteria", "Verrucomicrobia", "Lentisphaerae", "Fusobacteria", "Elusimicrobia", "Spirochaetes", "Synergistetes", "Tenericutes"))
prevdf.d7$Phylum <- factor(prevdf.d7$Phylum, levels = c("Firmicutes", "Bacteroidetes", "Proteobacteria", "Actinobacteria", "Verrucomicrobia", "Lentisphaerae", "Fusobacteria", "Elusimicrobia", "Spirochaetes", "Synergistetes", "Tenericutes"))

#Prevalence plot - t1
prevdf1.d_9 <- subset(prevdf.d_9, Phylum %in% get_taxa_unique(ps1.d_9, "Phylum"))
p.prevdf1.d_9 <- ggplot(prevdf1.d_9, aes(TotalAbundance, Prevalence / nsamples(ps1.d_9),color=Family)) +
  # geom_hline(yintercept = 0.03, alpha = 0.5, linetype = 2) +
  geom_point(size = 1.5, alpha = 0.7) +
  scale_x_log10(limits = c(1,10000), breaks = c(100,10000)) +
  xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
  facet_grid(~Phylum) +
  theme(legend.position="None") +
  theme(axis.title.x = element_blank()) +
  ggtitle("Day -9")

#Prevalence plot - t2
prevdf1.d0 <- subset(prevdf.d0, Phylum %in% get_taxa_unique(ps1.d0, "Phylum"))
p.prevdf1.d0 <- ggplot(prevdf1.d0, aes(TotalAbundance, Prevalence / nsamples(ps1.d0),color=Family)) +
  # geom_hline(yintercept = 0.03, alpha = 0.5, linetype = 2) +
  geom_point(size = 1.5, alpha = 0.7) +
  scale_x_log10(limits = c(1,10000), breaks = c(100,10000)) +
  xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
  facet_grid(~Phylum) +
  theme(legend.position="None") +
  theme(axis.title.x = element_blank()) +
  ggtitle("Day 0")

#Prevalence plot - t3
prevdf1.d7 <- subset(prevdf.d7, Phylum %in% get_taxa_unique(ps1.d7, "Phylum"))
p.prevdf1.d7 <- ggplot(prevdf1.d7, aes(TotalAbundance, Prevalence / nsamples(ps1.d7),color=Family)) +
  # geom_hline(yintercept = 0.03, alpha = 0.5, linetype = 2) +
  geom_point(size = 1.5, alpha = 0.7) +
  scale_x_log10(limits = c(1,10000), breaks = c(100,10000)) +
  xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
  facet_grid(~Phylum) +
  theme(legend.position="None") +
  theme(axis.title.x = element_blank()) +
  ggtitle("Day 7")

grid.arrange(p.prevdf1.d_9, p.prevdf1.d0, p.prevdf1.d7, nrow = 3)

Interpretation: There are a number of low prevalent bacteria. There is a notable reduction in both Firmicutes and Bacteroidetes following antibiotic treatment. This is not subset by treatment (randomization_arm), but the overall reduction suggests the antibiotics reduced the prevalence of many taxa following treatment. There is a visible restoration by day 7.

Interpretation: Log10 transformation does tend to shift the data to normality, but it is not striking and may ‘overshift’ the data. Should consider using log10 transformed data for future analysis, but consider other transformations as needed.

Community composition plotting

# Community composition - Baseline and Week 24 timepoints combined
# Phyla level plots
# Melt to long format (for ggploting) 
# Prune out phyla below % in each sample

# Creat a Phyla level table for external analysis
# write.table(ps1_phylum, file = "./results/phylum.txt", sep="\t")

threshold = 0.03

ps1_phylum <- ps1 %>%
  tax_glom(taxrank = "Phylum") %>%                     # agglomerate at phylum level
  transform_sample_counts(function(x) {x/sum(x)} ) %>% # Transform to rel. abundance
  psmelt() %>%                                         # Melt to long format
  filter(Abundance > threshold) %>%                    # Filter out low abundance taxa
  arrange(Phylum)                                      # Sort data frame alphabetically by phylum

# Convert Sample No to a factor because R is weird sometime
ps1_phylum$patient_ID <- as.factor(ps1_phylum$patient_ID)

p.comm.bar <- ggplot(ps1_phylum, aes(x = patient_ID, y = Abundance, fill = Phylum)) + 
  geom_bar(stat = "identity", width = 0.9) +
  facet_wrap(randomization_arm~day, scales = "free_x") +
  ggtitle("Community Composition Over Time and Randomization Arm") +
  ylab("Relative Abundance") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_blank()) +
  theme(axis.title.x = element_blank()) +
  scale_fill_brewer(palette = "Dark2")
p.comm.bar

# Draw interactive plot
ggplotly(p.comm.bar)

Interpretation: Similar to the prevalence plots, there is a notable alteration in the Bacteroidetes and Firmicutes in antibiotic treated samples at Day 0. This is still noticeable at Day 7, but the Firmicutes are making somewhat of a comeback. Proteobacteria invade following antibiotic treatment as well, but are largely absent, particularly in the narrow spectrum antibiotic samples after 7 days. Overall, the restoration to the pre-antibiotic bacterial microbiome appears to be more prominent following narrow spectrum treatment than in broad spectrum treatment.

Alpha diversity - categorical analysis

my.comparisons.paired <- list(c("Pre", "0"), c("Pre", "7"), c("0","7"))

p.paired.rich <- ggpaired(ps1.rich, x = "day", y = "richness_0", outlier.shape = NA, id = "patient_ID") +
  geom_jitter(width = 0.2) +
  ylab("Richness") +
  theme(axis.title.x = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  facet_grid(~randomization_arm) +
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = "-9", hide.ns = TRUE) +
  theme(axis.text.x = element_blank())

p.paired.sd <- ggpaired(ps1.rich, x = "day", y = "diversities_shannon", outlier.shape = NA, id = "patient_ID") +
  geom_jitter(width = 0.2) +
  ylab("Shannon Diversity") +
  theme(axis.title.x = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  facet_grid(~randomization_arm) +
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = "-9", hide.ns = TRUE)

grid.arrange(p.paired.rich, p.paired.sd, nrow = 2)

Interpretation: Richness (the number of taxa) and Shannon diversity (relative abundance and evenness) are similar in all untreated sampels over time. Both narrow and broad spectrum antibiotic treatment decrease alpha diversity over time, with some recovery occurring in both treatement protocols between Day 0 and 7.

Alpha diversity and day 7 rotavirus IgA boost

IgA boost defined as a greater than or equal to 2 fold change between pre-rotavirus IgA levels and Day 7 IgA levels.

p.rich.rota_boost <- ggboxplot(ps1.rich, x = "d7_rota_boost_updated", y = "richness_0", outlier.shape = NA) +
  geom_jitter(width = 0.2) +
  ylim(0,250) +
  ylab("Richness") +
  theme(axis.title.x = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  stat_compare_means(label = "p.signif", method = "t.test") +
  facet_grid(randomization_arm~day) +
  ggtitle("Rota-IgA Boost") +
  theme(axis.text.x = element_blank())

p.sd.rota_boost <- ggboxplot(ps1.rich, x = "d7_rota_boost_updated", y = "diversities_shannon", outlier.shape = NA) +
  geom_jitter(width = 0.2) +
  ylim(0,5) +
  ylab("Shannon diversity") +
  theme(axis.title.x = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  stat_compare_means(label = "p.signif", method = "t.test") +
  facet_grid(randomization_arm~day)

grid.arrange(p.rich.rota_boost, p.sd.rota_boost, nrow = 2)

Interpretation: There are signficiantly more types of bacteria (richness) at day 7 in both narrow spectrum and broad spectrum antibiotic patients. There is increased bacterial diversity in day -9 pre-treatment samples in the patients that go on to not receive antibiotics, but the levels are normalized over the course of the study.

Alpha diversity and shedding

Shedding was defined as having one or more stool samples per patient positive for rotavirus shedding.

p.rich.shedding <- ggboxplot(ps1.rich, x = "Shedding", y = "richness_0", outlier.shape = NA) +
  geom_jitter(width = 0.2) +
  ylim(0,250) +
  ylab("Richness") +
  theme(axis.title.x = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  stat_compare_means(label = "p.signif", method = "t.test") +
  facet_grid(randomization_arm~day) +
  ggtitle("Shedding") +
  theme(axis.text.x = element_blank())

p.sd.shedding <- ggboxplot(ps1.rich, x = "Shedding", y = "diversities_shannon", outlier.shape = NA) +
  geom_jitter(width = 0.2) +
  ylim(0,5) +
  ylab("Shannon diversity") +
  theme(axis.title.x = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  stat_compare_means(label = "p.signif", method = "t.test") +
  facet_grid(randomization_arm~day)

grid.arrange(p.rich.shedding, p.sd.shedding, nrow = 2)

Interpretation: Both richenss and diversity were increased in day 0 samples of patients treated with narrow spectrum antibiotics.

ANCOVA analysis

# Test for age dependency on alpha diversity
# Age per randomization_arm
ancova(richness_0~Age*randomization_arm, data = ps1.rich)
## Analysis of Variance Table
## 
## Response: richness_0
##                       Df Sum Sq Mean Sq F value   Pr(>F)   
## Age                    1     87    87.3  0.0309 0.860816   
## randomization_arm      2  28660 14329.9  5.0724 0.008092 **
## Age:randomization_arm  2   3549  1774.6  0.6282 0.535799   
## Residuals             94 265557  2825.1                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ancova(diversities_shannon~Age*randomization_arm, data = ps1.rich)
## Analysis of Variance Table
## 
## Response: diversities_shannon
##                       Df Sum Sq Mean Sq F value    Pr(>F)    
## Age                    1  0.809  0.8093  1.1815 0.2798223    
## randomization_arm      2 13.361  6.6806  9.7530 0.0001416 ***
## Age:randomization_arm  2  0.296  0.1480  0.2160 0.8061022    
## Residuals             94 64.388  0.6850                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# BMI
ancova(richness_0~BMIcalc*randomization_arm, data = ps1.rich)
## Analysis of Variance Table
## 
## Response: richness_0
##                           Df Sum Sq Mean Sq F value  Pr(>F)  
## BMIcalc                    1   2372  2372.3  0.8461 0.36000  
## randomization_arm          2  27040 13520.0  4.8222 0.01015 *
## BMIcalc:randomization_arm  2   4891  2445.3  0.8722 0.42140  
## Residuals                 94 263550  2803.7                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ancova(diversities_shannon~BMIcalc*randomization_arm, data = ps1.rich)
## Analysis of Variance Table
## 
## Response: diversities_shannon
##                           Df Sum Sq Mean Sq F value    Pr(>F)    
## BMIcalc                    1  0.372  0.3717  0.5545 0.4583477    
## randomization_arm          2 13.049  6.5243  9.7334 0.0001439 ***
## BMIcalc:randomization_arm  2  2.425  1.2126  1.8090 0.1694713    
## Residuals                 94 63.009  0.6703                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Alpha diversity is dependent on randomization arm, but not BMI

Interpretation: Alpha diversity is not dependent on age or BMI.

## All sample - Treatment
# Richness
p.gam.rich.treatment <- ggplot(ps1.rich, aes(x = day, y = richness_0, color=d7_rota_boost_updated)) +
  geom_point(alpha = 0.6, size = 3) +
  stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 3)) +
  ggtitle("General Additive Model (GAM) - Treatment") +
  facet_grid(~randomization_arm) +
  ylab("Richness")

# Shannon Diversity
p.gam.sd.treatment <- ggplot(ps1.rich, aes(x = day, y = richness_0, color=d7_rota_boost_updated)) +
  geom_point(alpha = 0.6, size = 3) +
  stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 3)) +
  ggtitle("General Additive Model (GAM) - Treatment") +
  facet_grid(~randomization_arm) +
  ylab("Shannon Diversity")

grid.arrange(p.gam.rich.treatment, p.gam.sd.treatment, nrow = 2)

## All sample - Treatment
# Richness
p.gam.rich.shedding <- ggplot(ps1.rich, aes(x = day, y = richness_0, color = Shedding)) +
  geom_point(alpha = 0.6, size = 3) +
  stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 3)) +
  ggtitle("General Additive Model (GAM) - Treatment") +
  facet_grid(~randomization_arm) +
  ylab("Richness")

# Shannon Diversity
p.gam.sd.shedding <- ggplot(ps1.rich, aes(x = day, y = richness_0, color=Shedding)) +
  geom_point(alpha = 0.6, size = 3) +
  stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 3)) +
  ggtitle("General Additive Model (GAM) - Treatment") +
  facet_grid(~randomization_arm) +
  ylab("Shannon Diversity")

grid.arrange(p.gam.rich.shedding, p.gam.sd.shedding, nrow = 2)

Beta diversity analysis

# Ordination

# weighted unifrac
ord.ps1.wuni <- ordinate(ps1, method = "PCoA", distance = "wunifrac")
ord.ps1.wuni.log <- ordinate(ps1.log, method = "PCoA", distance = "wunifrac")

# unifrac
ord.ps1.uni <- ordinate(ps1, method = "PCoA", distance = "unifrac")

# Bray-curtis
ord.ps1.bray <- ordinate(ps1, method = "PCoA", distance = "bray")

# PCoA plot
p.ord.wuni <- plot_ordination(ps1, ord.ps1.wuni, color = "randomization_arm") +
  geom_point(size=3, alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  geom_point(colour = "grey50", size = 1.5) +
  facet_grid(~ day) +
  stat_ellipse(geom = "polygon", alpha = 1/15, aes(fill = randomization_arm)) +
  labs(color = "randomization_arm") +
  ggtitle("wUniFrac")

p.ord.wuni.log <- plot_ordination(ps1.log, ord.ps1.wuni.log, color = "randomization_arm") +
  geom_point(size=3, alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  geom_point(colour = "grey50", size = 1.5) +
  facet_grid(~ day) +
  stat_ellipse(geom = "polygon", alpha = 1/15, aes(fill = randomization_arm)) +
  labs(color = "randomization_arm") +
  ggtitle("wUniFrac: log transformed")

p.ord.uni <- plot_ordination(ps1, ord.ps1.uni, color = "randomization_arm") +
  geom_point(size=3, alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  geom_point(colour = "grey50", size = 1.5) +
  facet_grid(~ day) +
  stat_ellipse(geom = "polygon", alpha = 1/15, aes(fill = randomization_arm)) +
  labs(color = "randomization_arm") +
  ggtitle("UniFrac")

p.ord.bray <- plot_ordination(ps1, ord.ps1.bray, color = "randomization_arm") +
  geom_point(size=3, alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  geom_point(colour = "grey50", size = 1.5) +
  facet_grid(~ day) +
  stat_ellipse(geom = "polygon", alpha = 1/15, aes(fill = randomization_arm)) +
  labs(color = "randomization_arm") +
  ggtitle("Bray-Curtis")

grid.arrange(p.ord.wuni, p.ord.wuni.log, p.ord.uni, p.ord.bray, nrow = 4)

Interpretation: Selecting wUniFrac as it 1) tends to agree with alpha diversity observations and 2) explains a high-percentage of variance on axis 1. However, there are likely details to be resolved from each measure.

Community composition differences beetween groups is primarily explained by treatment (randomization arm). As with alpha diversity, the extent of difference between groups is greatest at day 0 and diminished by day 7. However, the differences between each sample at day 7 are larger than at day 7 (more variant across Axis 1 and 2) suggesting intra-individual spread in bacterial communities.

Examine beta diversity grouping based on rotavirus boosting or rotavirus shedding in the stool.

p.ord.boost <- plot_ordination(ps1, ord.ps1.bray, color = "d7_rota_boost_updated") +
  geom_point(size=3, alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  geom_point(colour = "grey50", size = 1.5) +
  facet_grid(randomization_arm ~ day) +
  stat_ellipse(geom = "polygon", alpha = 1/15, aes(fill = d7_rota_boost_updated)) +
  ggtitle("wUniFrac Rotavirus IgA Boost")

p.ord.shedding <- plot_ordination(ps1, ord.ps1.bray, color = "Shedding") +
  geom_point(size=3, alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  geom_point(colour = "grey50", size = 1.5) +
  facet_grid(randomization_arm ~ day) +
  stat_ellipse(geom = "polygon", alpha = 1/15, aes(fill = Shedding)) +
  ggtitle("wUniFrac Rotavirus Shedding")

grid.arrange(p.ord.boost, p.ord.shedding, nrow = 2)

Interpretation: Visual inspection makes it appear that bacterial community structure is different between samples from patients that boosted and those that did not boost at Day 7. There is poor sample distribution of boosters and shedders to analyze no antibiotic controls and too few boosted samples to analyze broad spectrum samples.

Premanova significance testing (ADONIS)

# Set a random seed so that exact results can be reproduced
set.seed(1000)

# Function to run adonis test on a physeq object and a variable from metadata
doadonis <- function(physeq, category) {
  bdist <- phyloseq::distance(physeq, "wunifrac")
  col <- as(sample_data(physeq), "data.frame")[ ,category]
  
  # Adonis test
  adonis.bdist <- adonis(bdist ~ col)
  print("Adonis results:")
  print(adonis.bdist)
  
  # Homogeneity of dispersion test
  betatax = betadisper(bdist,col)
  p = permutest(betatax)
  print("Betadisper results:")
  print(p$tab)
}

# Two factor adonis
doadonis2 <- function(physeq, category, category2) {
  bdist <- phyloseq::distance(physeq, "wunifrac")
  col <- as(sample_data(physeq), "data.frame")[ ,category]
  col2 <- as(sample_data(physeq), "data.frame")[ ,category2]
  
  # Adonis test
  adonis.bdist <- adonis(bdist ~ col * col2)
  print("Adonis results:")
  print(adonis.bdist)

}

# Runs Permanov and Homogeneity of dispersion test
# See ?betadisper for more info

# Randomization arm at each time point
doadonis2(ps1, "day", "randomization_arm") # day = **, randomization_arm = ***, interaction = n.s.
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col * col2) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## col        1    0.5808 0.58085  3.9030 0.03398  0.005 ** 
## col2       2    2.0554 1.02770  6.9055 0.12025  0.001 ***
## col:col2   2    0.4665 0.23323  1.5672 0.02729  0.106    
## Residuals 94   13.9894 0.14882         0.81847           
## Total     99   17.0921                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
doadonis(ps1.t1, "randomization_arm") # n.s.
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)
## col        2   0.06800 0.033998  1.1857 0.10146  0.316
## Residuals 21   0.60215 0.028674         0.89854       
## Total     23   0.67014                  1.00000       
## [1] "Betadisper results:"
##           Df     Sum Sq     Mean Sq        F N.Perm Pr(>F)
## Groups     2 0.02441588 0.012207942 2.903593    999  0.083
## Residuals 21 0.08829296 0.004204427       NA     NA     NA
doadonis(ps1.t2, "randomization_arm") # ***
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## col        2   0.78594 0.39297  16.867 0.52112  0.001 ***
## Residuals 31   0.72223 0.02330         0.47888           
## Total     33   1.50817                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
##           Df      Sum Sq      Mean Sq         F N.Perm Pr(>F)
## Groups     2 0.001956101 0.0009780507 0.1497809    999  0.853
## Residuals 31 0.202426115 0.0065298747        NA     NA     NA
doadonis(ps1.t3, "randomization_arm") # **
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## col        2    0.6298 0.31489  2.5571 0.11593  0.019 *
## Residuals 39    4.8025 0.12314         0.88407         
## Total     41    5.4323                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
##           Df     Sum Sq    Mean Sq      F N.Perm Pr(>F)
## Groups     2 0.06811062 0.03405531 2.5677    999  0.088
## Residuals 39 0.51725553 0.01326296     NA     NA     NA
# Randomization arm over time
doadonis(ps1.con, "day") # n.s.
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)
## col        1   0.06192 0.061919  1.0496 0.03176  0.366
## Residuals 32   1.88783 0.058995         0.96824       
## Total     33   1.94974                  1.00000       
## [1] "Betadisper results:"
##           Df     Sum Sq     Mean Sq        F N.Perm Pr(>F)
## Groups     2 0.02073363 0.010366814 1.086088    999  0.345
## Residuals 31 0.29589795 0.009545095       NA     NA     NA
doadonis(ps1.narrow, "day") # n.s. (p < 0.1) # likley not-significant due to Day -9 and Day 7 similarity
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)
## col        1    0.0931 0.093058 0.90936 0.02605  0.427
## Residuals 34    3.4793 0.102333         0.97395       
## Total     35    3.5724                  1.00000       
## [1] "Betadisper results:"
##           Df    Sum Sq     Mean Sq        F N.Perm Pr(>F)
## Groups     2 0.1248340 0.062417018 9.693034    999  0.001
## Residuals 33 0.2124992 0.006439368       NA     NA     NA
doadonis(ps1.broad, "day") # ***
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## col        1    0.2855 0.28545  2.0451 0.06807  0.087 .
## Residuals 28    3.9082 0.13958         0.93193         
## Total     29    4.1936                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
##           Df    Sum Sq    Mean Sq        F N.Perm Pr(>F)
## Groups     2 0.4770339 0.23851694 20.64645    999  0.001
## Residuals 27 0.3119159 0.01155244       NA     NA     NA
# Interactions between day and rotavirus boost in separate treatment groups
doadonis2(ps1.narrow, "day", "d7_rota_boost_updated") # n.s.
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col * col2) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## col        1    0.1429 0.14286 0.83695 0.02379  0.500
## col2       1    0.2141 0.21414 1.25454 0.03566  0.253
## col:col2   1    0.1854 0.18545 1.08646 0.03088  0.318
## Residuals 32    5.4621 0.17069         0.90966       
## Total     35    6.0045                 1.00000
# Rotavirus sheeding in each radnomization arm over time
doadonis2(ps1.narrow, "day", "Shedding") # n.s.
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col * col2) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)
## col        1    0.1110 0.110998 0.68401 0.02035  0.603
## col2       1    0.0678 0.067804 0.41783 0.01243  0.814
## col:col2   1    0.0820 0.081964 0.50509 0.01503  0.752
## Residuals 32    5.1928 0.162276         0.95218       
## Total     35    5.4536                  1.00000
doadonis2(ps1.broad, "day", "Shedding") # *** for day, n.s. for shedding (samples are different over time, but not based on rotavirus shedding)
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col * col2) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## col        1   0.35058 0.35058  4.2481 0.13581  0.003 **
## col2       1   0.03771 0.03771  0.4569 0.01461  0.834   
## col:col2   1   0.04753 0.04753  0.5760 0.01841  0.760   
## Residuals 26   2.14568 0.08253         0.83117          
## Total     29   2.58150                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Interaction between treatment and rotavirus vaccine boost at each time point
doadonis2(ps1.t1, "randomization_arm", "d7_rota_boost_updated") # n.s.
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col * col2) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)
## col        2   0.17375 0.086873 1.20147 0.10114  0.295
## col2       1   0.06167 0.061674 0.85295 0.03590  0.422
## col:col2   1   0.10859 0.108588 1.50178 0.06321  0.184
## Residuals 19   1.37382 0.072306         0.79974       
## Total     23   1.71783                  1.00000
doadonis2(ps1.t2, "randomization_arm", "d7_rota_boost_updated") # n.s. (samples are differnet by randomization arm only)
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col * col2) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## col        2    3.5305 1.76524 16.9073 0.52575  0.001 ***
## col2       1    0.0562 0.05624  0.5387 0.00838  0.720    
## col:col2   1    0.1006 0.10061  0.9637 0.01498  0.447    
## Residuals 29    3.0278 0.10441         0.45089           
## Total     33    6.7152                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
doadonis2(ps1.t3, "randomization_arm", "d7_rota_boost_updated") # sasmples are differnt by randomization arm (**) and by rotavirus boost status, but there is no interaction
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col * col2) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## col        2   0.26196 0.13098  2.4310 0.09178  0.037 * 
## col2       1   0.58868 0.58868 10.9262 0.20626  0.002 **
## col:col2   2   0.06386 0.03193  0.5926 0.02238  0.589   
## Residuals 36   1.93960 0.05388         0.67958          
## Total     41   2.85409                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Double check day 7 for narrow as it looks separated in the plot
ps1.t3.narrow <- subset_samples(ps1.t3, randomization_arm == "Narrow spectrum antibiotics")
sample_data(ps1.t3.narrow)$randomization_arm
##  [1] Narrow spectrum antibiotics Narrow spectrum antibiotics
##  [3] Narrow spectrum antibiotics Narrow spectrum antibiotics
##  [5] Narrow spectrum antibiotics Narrow spectrum antibiotics
##  [7] Narrow spectrum antibiotics Narrow spectrum antibiotics
##  [9] Narrow spectrum antibiotics Narrow spectrum antibiotics
## [11] Narrow spectrum antibiotics Narrow spectrum antibiotics
## [13] Narrow spectrum antibiotics Narrow spectrum antibiotics
## [15] Narrow spectrum antibiotics
## Levels: Narrow spectrum antibiotics
sample_data(ps1.t3.narrow)$day
##  [1] 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
doadonis(ps1.t3.narrow, "d7_rota_boost_updated") # 0.08. < 0.1 is near significance
## [1] "Adonis results:"
## 
## Call:
## adonis(formula = bdist ~ col) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## col        1    0.3965 0.39650  2.6446 0.16904  0.087 .
## Residuals 13    1.9491 0.14993         0.83096         
## Total     14    2.3456                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
##           Df     Sum Sq    Mean Sq        F N.Perm Pr(>F)
## Groups     1 0.01071119 0.01071119 0.307049    999  0.571
## Residuals 13 0.45349577 0.03488429       NA     NA     NA

Interpretation: